/////////////////////////////////////////////////////////////////////
//
// \author Jeongnim Kim  jnkim@ncsa.uiuc.edu
//
// Measuring abstraction penalty using TinyVector and Vector based on PETE
// src/PETE      (copy of PETE library http://www.acl.lanl.gov)
// src/OhmmsPETE (Ohmms implementation of a fixed vector and an array)
//////////////////////////////////////////////////////////////////////
#include <iostream>
#include <vector>
#include <iomanip>
using namespace std;

#include "Pooma/Arrays.h"
#include "Utilities/Clock.h"
 
#if defined(ADD_)
#define addmult addmult_
#endif

extern "C" { // fortran call
  void addmult(double* a,const double* b,const double* c, const double *d,
	       int* ntot);
}

inline void 
caddmult(double* restrict a, const double* restrict b, 
	 const double* restrict c, 
	 const double* restrict d, int ntot) 
{
  for(int i=0; i<ntot; i++) {
     a[i] = b[i] + 1.3*c[i] + 1.7*d[i];
  } 
}


int main(int argc, char **argv) {

   // Initialize Pooma.
  Pooma::initialize(argc, argv);

  typedef Vector<3> Pos_t;
  typedef Array<1, Vector<3>, MultiPatch<UniformTag,Brick>  > PosArray_t;

  int n=10;
  int i=0;
  while(i<argc) {
    if(!strcmp(argv[i],"-n")) {
      n = atoi(argv[++i]);
      break;
    }
    i++;
  }
  

  cout << endl << "Starting vector operations:" << n << " " << OHMMS_DIM
       <<"-dimensional elements" << endl;

  Loc<1> patchSize(10);
  UniformGridLayout<1> layout(Interval<1>(n), patchSize, ReplicatedTag());

  Pooma::Clock myclock;
  double tstart, tend;
  vector<double> timingdata(10);

  tstart = myclock.value();
  PosArray_t A(layout), B(layout), C(layout), D(layout);
  tend = myclock.value();
  timingdata[0] = tend-tstart;

  tstart = myclock.value();
  A = 0.0;
  for(int i=0; i<n; i++) {
    B(i) = Pos_t(2*rand()-1, 2*rand()-1, 2*rand()-1);
  }
  for(int i=0; i<n; i++) {
    C(i) = Pos_t(2*rand()-1, 2*rand()-1, 2*rand()-1);
  }
  D = 1.1*B;

  tend = myclock.value();
  timingdata[1] = tend-tstart;

  tstart = myclock.value();
  A = B + 1.3*C + 1.7*D;
  tend = myclock.value();
  timingdata[2] = tend-tstart;


  tstart = myclock.value();
  for(int i=0; i<n; i++)
    A(i) = B(i) + 1.3*C(i) + 1.7*D(i);
  tend = myclock.value();
  timingdata[3] = tend-tstart;

 
  if(n<=100) {
    for(i=0; i<n; i++)   {
      cout << A(i) << " " << B(i) << " " << C(i) <<  " " << D(i) << endl;
    }
  }

  typedef Array<1,double> ValArray_t;
  ValArray_t a(OHMMS_DIM*n), b(OHMMS_DIM*n), c(OHMMS_DIM*n), d(OHMMS_DIM*n);
  int ii=0;
  for(int i=0; i<n; i++) {
    for(int j=0; j<OHMMS_DIM; j++)
      b(ii++) = B(i)(j);
  }

  ii=0;
  for(int i=0; i<n; i++) {
    for(int j=0; j<OHMMS_DIM; j++)
      c(ii++) = C(i)(j);
  }

  ii=0;
  for(int i=0; i<n; i++) {
    for(int j=0; j<OHMMS_DIM; j++)
      d(ii++) = D(i)(j);
  }

  tstart = myclock.value();
  a = b + 1.3*c + 1.7*d;
  tend = myclock.value();
  timingdata[4] = tend-tstart;


  int ntot = OHMMS_DIM*n;
  tstart = myclock.value();
  for(int i=0; i<ntot; i++) a(i) = b(i) + 1.3*c(i) + 1.7*d(i);
  tend = myclock.value();
  timingdata[5] = tend-tstart;


  tstart = myclock.value();
  caddmult(&(a(0)), &(b(0)), &(c(0)), &(d(0)), ntot);
  tend = myclock.value();
  timingdata[6] = tend-tstart;


  tstart = myclock.value();
  addmult(&(a(0)), &(b(0)), &(c(0)), &(d(0)), &ntot);
  tend = myclock.value();
  timingdata[7] = tend-tstart;

  cout.setf(ios::fixed, ios::floatfield);
  cout.setf(ios::right, ios::adjustfield);

  cout << "Vector<TinyVector<3> > Allocation " << setw(10) << timingdata[0] << endl;
  cout << "Vector<TinyVector<3> > Assignment " << setw(10) << timingdata[1] << endl;

  cout << endl
       << "Operations A = B + 1.3*C + 1.7* D  "
	<< setw(10) << "CPU secs" << setw(10) << "Penalty"
        << " C-penalty" << endl;
  cout << "AbstractionTest0           " 
       << setw(10) << timingdata[2] 
       << setw(10) << timingdata[2]/timingdata[5] 
       << setw(10) << timingdata[2]/timingdata[7]<< endl;

  cout << "AbstractionTest1           "
       << setw(10) << timingdata[3] 
       << setw(10) << timingdata[3]/timingdata[5] 
       << setw(10) << timingdata[3]/timingdata[7]<< endl;

  cout << "AbstractionTest2           " 
       << setw(10) << timingdata[4] 
       << setw(10) << timingdata[4]/timingdata[5] 
       << setw(10) << timingdata[4]/timingdata[7]<< endl;

  cout << "AbstractionTest3           " 
       << setw(10) << timingdata[5] 
       << setw(10) << timingdata[5]/timingdata[5] 
       << setw(10) << timingdata[5]/timingdata[7]<< endl;

  cout << "C inline function restrict " 
       << setw(10) << timingdata[6] 
       << setw(20) << timingdata[6]/timingdata[7]<< endl;
  cout << "fortran subroutine         " 
       << setw(10) << timingdata[7] 
       << setw(20) << timingdata[7]/timingdata[7]<< endl;

  cout << endl;
}

/***************************************************************************
 * $RCSfile: poomavector.cpp,v $   $Author: jnkim $
 * $Revision: 1.1.1.1 $   $Date: 2004/07/23 21:04:02 $
 * $Id: poomavector.cpp,v 1.1.1.1 2004/07/23 21:04:02 jnkim Exp $
 ************************************************************************/
